Critical Scaling in Linear Response of Prictionless Granular Packings near Jamming 
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We study the origin of the scaling behavior in frictionless granular media above the jamming 
transition by analyzing their linear response. The response to local forcing is non-self-averaging 
and fluctuates over a length scale that diverges at the jamming transition. The response to global 
forcing becomes increasingly non-afhne near the jamming transition. This is due to the proximity 
of floppy modes, the influence of which we characterize by the local linear response. We show that 
the local response also governs the anomalous scaling of elastic constants and contact number. 

PACS numbers: 45.70.-n, 46.65-fg, 83.80.Fg, 05.40.-a 



The general picture of jamming which was advanced 
for systems U y, S H that form a shear-resistent solid 
phase at high densities, is bringing a new perspective 
to the deformations of granular and disordered media. 
A good model for studying such media arepackings of 
polydisperse weakly compressible spheres |3, 14- If we 
measure pressure in units of the elastic constants and 
characteristic radius of the balls (as we will do below), 
the relevant limit for granulates is the small- deformation 
or, equivalently, the small-pressure limit in the absence 
of thermal fluctuations. This limit is also relevant for 
weakly compressed emulsions |5j- We will focus on the 
case of frictionless, deformable spherical particles, and 
introduce a simple, experimentally accessible and local 
measure to characterize the nature of their deformations 

Deformable particles form a stiff jammed phase when 
the pressure becomes larger than zero. At the zero pres- 
sure jamming point 'J', packings form a "marginal solid" 
and are isostatic, i.e, the average number of contacts per 
particle, z, reaches the minimum zPq = 2(i, needed for 
a frictionless packing to remain stable in d dimensions. 
When the point J is approached by decreasing the pres- 
sure, several surprising scaling relations emerge: the ex- 
cess contact number Az = z— zP^ scales as yS, with S the 
typical dimensionless compression of the particles, while 
the ratio G/K of the shear modulus G to the compres- 
sion modulus K scales as Az. In addition, a diverging 
time scale w* ~ Az has been identified in the density of 
states of vibrational modes. The jamming point J thus 
exhibits features of a critical point |^, y, |^ . 

Since packings at the jamming point are marginal, ev- 
ery additional broken contact generates a global zero- 
energy displacement mode, a so-called floppy mode [3,|3j 
M>M- Wyart and coworkers |alEll3 have shown that the 
scaling near J is related to those floppy modes, by creat- 
ing trial modes for the deformations of weakly jammed 
solids. These modes are based on the floppy modes that 
would occur when along the faces of cubes of linear size 
£* ~ 1/Az bonds would be cut. Even though for jammed 
systems, truly floppy modes never occur, their proximity 
thus governs the scaling just above the jamming point. 
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FIG. 1: (color) Force response networks for a point loading 
with pressure as indicated. Blue (red) lines indicate positive 
(negative) changes in contact force, the thickness indicating 
the amount. The particles themselves are not drawn. 



In this Letter we uncover that this proximity of floppy 
modes causes an increasingly non-affine response when 
approaching point J, and that this response is intimately 
related to the (anomalous) scalings of the shear modulus 
G, the excess contact number Az and the length scale 
£* . We numerically study the linear, quasistatic response 
of systems near the jamming transition. The resp onse of 
granular media has been widely-studied [III ll^ U^ llJi 
llSlllq . but not, we believe, systematically as a function 
of the distance to the jamming point J. Nor does it seem 
to have been fully appreciated that the scaling behavior 
can essentially be captured within linear response. 

We represent the linear response by relative displace- 
ments and changes in contact forces, and find significant 
changes with the distance to point J. (i) Fig. 1 illustrates 
that the response to the loading of a single grain be- 
comes increasingly disordered over an increasingly large 
scale when the jamming transition is approached — this 
leads to a direct observation of the diverging length scale 
£* ~ 1/Az, shown below. We will show that such a lo- 
cal force response is not self-averaging, even though it 
is smooth upon ensemble-averaging and then quantita- 
tively agrees with continuum elastic behavior, (ii) The 
response to a uniformly applied compression or shear also 
varies with the distance to jamming. We introduce the 
distribution P{a) of angles a between the bonds and the 



local deformations as an indicator of the non-afRne na- 
ture of the response. Near J, P{a) becomes strongly 
peaked around a — tt/2, with the width and height of 
the peak scaling with the distance to the jamming point. 
Grains then predominantly slide past each other, which 
signals an increasingly non-afhne response of the mate- 
rial caused by the proximity of floppy modes, for which 
P{a) = d{a — tt/2). (in) The component of the relative 
displacements perpendicular to the bond vector diverges 
upon approaching the jamming point, (iv) Finally, the 
Az ~ V6 scaling j3j is identified to originate precisely 
from this increasingly "sliding" response. 

Hence a simple picture emerges: the influence of floppy 
modes can be quantifled by the local linear response of 
the material, which becomes increasingly non-afflne near 
jamming, in turn causing anomalous scaling. 

Linear response — The response of a jammed granular 
medium to external loads has been studied mostly by full 
scale molecular dynamics [13, llJI • We calculate it here 
in linear order from expansion of energy to 2"*^ order 
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Here the sum is over all contacts, Ui is the displacement 
of particle i, and Uy = Uj — Ui the relative displacement 
of grains i and j, with components u\\^ij and u^^ij parallel 



and perpendicular to the bond vector r^ 
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notes the stiffness of the contact and /^ its initial force. 
The second term proportional to u\ - ■ is due to the trans- 
verse motion when the bonds are pre-stressed {fij 7^ 0). 
For contact interactions which increase as a power of the 
overlap fij ~ S'^^, the factor fij/hjrij = Sij/{l3rij) is of 
order of the dimcnsionless compression d — Sij/rij, which 
is small and which vanishes at the jamming point. 

We study 2D packings of N frictionless Hertzian 
3/2 
spheres for which fij ^ 6,^' , where 6ij is the overlap 

between neighboring particles. The conflning pressure 



ranges from p = 10 to p 



lo- 



in units of the ef- 



fective Young modulus of the constituent particles. See 
Ref. [13 for details. For each packing, the expansion (1) 
yields the dynamical matrix M. Instead of studying the 
vibrational dynamics [a, |3, |l3, llSj we obtain here the 
quasistatic response to external forces f^^^ |20j by solv- 
ing the linear equation Mij^af3Uj,i3 = ff^, for the Uj and, 
through the force law, the forces fij. Here i,j label the 
particles and a, /? the coordinate axes. 

Elastic Moduli — We have calculated the elastic mod- 
uli from the linear response by applying an overall com- 
pression or shear and by point-loading a single particle, 
for packings with N = 10"^ and N = 10"* respectively. 
The resulting force fields are translated into local stress 
fields llfl] which are then ensemble averaged. From fits of 
the point response [21J , we determine both K and G and 
compare these to the values obtained from the response 
to global shear and compression. Fig. 2a shows that these 
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FIG. 2; (a) Scaling of bulk/shear modulus with pressure, ob- 
tained from point response (squares), and global compres- 
sion/shear (diamonds). The fitted exponents are 0.38 ± 0.03 
for the bulk and 0.70 ± 0.08 for the shear modulus, (b) Dis- 
tribution function of the scaled transverse response for shear 
(see text), (c-d) Distribution of the relative displacement an- 
gle a for (c) shear and (d) compression of packings for a range 
of pressures. Insets: scaling of the width of the peak ~ p^ . 



two methods agree very well quantitatively, and that the 
elastic moduli scale with pressure as i^ ~ P^^'^, G ^ p"^^^, 



in agreement with earlier results [3, 122 

Non-affinity — The typical bond stiffness kij is pro- 
portional to p^/^ for Hertzian contacts. Hence, a simple 
estimate for the elastic moduli scaling as p^^"^ follows un- 
der the affinity assumption that the bond deformations 
are of order of the applied deformation. This estimate 
fails for the shear modulus G, which vanishes faster than 
K when approaching the jamming point: G/K ^ Az. 
This has been thought to be caused by strongly non-affine 
behavior of the system under shear J3| and the proximity 
of the floppy modes |[ly| . We will elucidate now the cause 
of the scaling and the influence of the floppy modes via 
the local deformations U||_y and u_\_,ij. 

As the eigenmodes or snapshots of the response look 
very disordered [ifl ll3) llM HOl i it has turned out to be 
difficult to find a simple measure to characterize the non- 
affinity and the overall floppy mode character. We show 
now that proximity of the floppy modes can clearly be 
identified in the distribution P{a) of the local angles 
aij = &ta,n{u\\ij / u^^ij) . In a disordered, isotropic sys- 
tem, one expects P{a) = 5{a — tt) for a purely homoge- 
nous compression, and P{a) = I/tt for a purely affine 
shear. In contrast, for a floppy mode, P{a) = 5{a — tt/2). 
This is because in floppy modes the relative angles be- 
tween particles change while the relative distances r^j 
remain unchanged, as if all the bonds are replaced by 
incompressible sticks Q- Hence, for a true floppy mode 
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As the pressure is lowered, P{a) and P{u±) evidence 
that the local deformations evolve from near-affine to ex- 
tremely non-affine, floppy-mode like (Fig. 2b-d). Indeed, 
for a sheared system, P{a) evolves from a flat distribu- 
tion at large pressures to a sharply peaked distribution 
for lower pressures (Fig. 2c). This peak is located around 
7r/2 and its weight approaches 1 — locally the response 
becomes more and more transverse for p —^ and P{a) 
approaches that of a floppy mode. For a compressed sys- 
tem at large pressures, P{a) has the "afflne" peak around 
a = TT, while for lower pressures -P(a) again develops a 
sharp peak around a — Tr/2 (Fig. 2d). 

Even though the response is far from afflne for both 
compression and shear, the afflne prediction for K holds 
true while it fails for G. The reason is that for compres- 
sion, only a flnite fraction of the displacements is essen- 
tially transverse and P{a) remains non-zero away from 
the peak at tt/2. Since according to the energy expres- 
sion (^ the compression of bonds given by u^^ij gives the 
dominant contribution to AE, this is consistent with the 
fact that the compression modulus scales with the bond 
stiffness k: K ^ k ^ p^^"^. For shear deformations, how- 
ever, fewer and fewer bonds contribute to leading order 
to the energy, and the weight outside the peak vanishes 
as Az ~ p^/^, consistent with the scaling G/K ^ Az. 

Scaling of P{a) and P{u±) — We can understand the 
development of the peak in P{a) from the balance of 
terms in the energy expansion (^. Focussing on typical 
values, AE ^ fc(u| — Su^). Since k ^ p^/^ and 5 ^ 



Az^ ~ p^l'^ ^ balancing the terms we flnd that 






nl/3 



(2) 



SO that for small p, P[a) develops a peak around a = 71/2, 
and the width of this peak should scale as p^/^. This is 
what we flnd — see the insets of Figs. 2c, d. 

How do the typical values U|| and ux scale when we im- 
pose a global shear or compression of order 7 on the sys- 
tem? Equating the elastic energy densities for compres- 
sion and shear, K^"^ and 67^, to the energy expansion, 
and knowing that the elastic moduli scale as G ~ p 
and K ^ p^^^, we can predict the scaling of uu and uj_: 
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These scaling predictions are well obeyed by our data for 
small p — in Fig. 2b we show this for shear deformations 
and P{u±). 

The fact that for flxed 7 the typical perpendicular re- 
sponse u± diverges upon approaching the jamming point, 
is connected to the disordered nature of the microscopic 
response already familiar from the randomly oriented 
swirl-type motions [llJ, UJ, lla] that characterize eigen- 
modes and responses. Of course, in a system of flnite 
size, uj_ can not diverge. If the cross-over is determined 



by the length-scale £* becoming of order the linear system 
size L, one expects a cross-over scaling uj_ = L^'^g{i* /L) 
with g{w) ^ w^'^ for w -^ 1 and g ~ const for w 3> 1. 
Note also that in the regime L/£* <C 1 the response is 
close to that of a floppy mode appearing at isostaticity, 
while our scaling results apply to the regime L/£* ^ 1. 

Az scaling — The non-afflne response also nicely ex- 
plains the microscopic origin of the anomalous Az ^ vS 
scaling under compression — theories assuming afflne de- 
formations give Az ~ d. Let us consider a small compres- 
sion of the packing with typical bond compression U||. 
This leads to an inflnitesimal change in contact number 
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Upon lowering the pressure, the global compression will 
excite distorted floppy modes, i.e., for many bonds u±^ij 
will be of order u\\/V6 (see Eq. (0) and Fig. 2d). More- 
over, the chance that in such almost perpendicular dis- 
placements one of the particles bumps into a nearby par- 
ticle with which it was not in contact yet, will be propor- 
tional to this motion, i.e. to u^^/y/S. Equating this to (O 
then yields Az ^ y/5. The picture which thus emerges is 
that the lower the pressure, the larger the u±^ij, and the 
larger the chances that new contacts are created. 

In earlier papers, it was noted [j, U3, l2^ that the 
Az ~ v5 scaling with compression was consistent with 
a square root divergent term in the correlation function 
g{r), if it was assumed that compression would be essen- 
tially afflne. As we have seen, however, distortions are 
not at all afflne near jamming. Our analysis turns this 
around: it suggests that the natural coordinates for the 
floppy mode-like distortions are the perpendicular dis- 
placements, not the radial ones, and that these generate 
the square root behavior of g{r) in the radial direction. 

Point response and diverging length scale — We flnally 
return to the point response. Fig. 3a illustrates that the 
ensemble average of such a response conforms to elastic- 
ity; the stress fields fit well and the fitted elastic constants 
agree with those obtained from bulk response (Fig. 2a). 
On the other hand, individual responses become very dis- 
ordered and suggest the occurrence of a large length scale 
£* when approaching point J (see Fig. 1). 

To extract this length scale, we characterize the re- 
sponse to an infinitesimal infiation of a single central 
grain, since the response to a local directional force, as 
shown in Fig. 1 and 3a, is highly anisotropic. We nor- 
malize the forces by fitting the radial stress to the elastic 
response [2g. We focus on the radial component of the 
change in contact force, dfr, calculate {dfr{r)) by averag- 
ing dfr over concentric rings, and study the RMS fluctu- 

ations h{r) ^ ^ {[dUr) -{df^))?)- 

The range over which these fluctuations are felt grows 
when approaching the jamming transition (Fig. 3b). 
When plotted as function of rAz, the data for h col- 
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FIG. 3: (a) Shear stress a^y 
and p = 10 
diameter [l( 

fitted to linear elasticity (dashed), (b-d) Spatial decay of 
the force fluctuations h. Squares, diamonds, crosses and plus 
signs represent increasing pressures, ranging from p = 4.6 ■ 
10~^ to p = 3.2 • 10~^, respectively, (b) h as function of r. 
(c) Scaling collapse when h is plotted as function of rAz. (d) 
Scaling collapse of relative fluctuations. 



lapses (Fig. 3c). To our knowledge, this is the first ev- 
idence for the existence of a length scale £* ^ 1/Az in 
local response measurements. Essentially the same char- 
acteristic length was shown by Wyart et al. to govern the 
vibrational density of states [8|,|9|. This scale is identified 
as the linear size i* of the largest domain which deforms 
freely by pushing on the bonds at its surface. Equat- 
ing the number of bonds on the surface (~ £'^~^) to the 
number of excess bonds in the bulk (~ Az £'') yields the 
maximum size of such domain i* ~ 1/Az [alSl- 

Both the average response {df^{r)) and the fluctua- 
tions h{r) decay as 1/r^ — the relative fluctuations do 
not decay far from the perturbed grain. The response 
is not self-averaging, and there is no finite correlation 
length of fluctuations. The asymptotic value of the rela- 
tive fluctuations h{r) / {df ^{r)) = h{r)r^ grows as 1/Az^ 
(Fig. 3d) — one has to coarse grain the response over 
increasingly more grains 0{{Az)~'^) to start to see con- 
vergence to average continuum-like stress response. 

Outlook — Our analysis of the (local) linear response 
substantiates and extends the concept that the jammed 
phase of weakly compressed frictionless particles is dom- 
inated by the proximity of floppy modes 0, H B llUl • We 
identified the increasingly non-afline response to give rise 
to the scaling Az ^ \/S and presented a direct observa- 
tion of the scale £* ^ l/Az introduced before 0,13 • The 
emerging scenario favors a microscopic, geometric inter- 
pretation of these scalings and has several implications 
that deserve further study: (i) What is the finite size 
scaling form of u±7 (ii) What happens for non-power- 
law contact interactions such as fij ~ exp(— (5^- ■" ') with 



/? > 1? Our analysis suggests a scaling Az ~ 5^^'^. 
(iii) What happens to the square root divergence in g(r) 
Pl I2J1 12^ when the packing algorithm does not allow for 
floppy mode- like rearrangements upon annealing, such as 
may occur for algorithms based on local rearrangements 
or packings of truly hard spheres? 

We thank K. Shundyak, J.H. Snoeijer, and M. Dep- 
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dation FOM and MvH support from NWO/VIDI. 



[1] 
[2] 
[3] 
[4] 
[5] 

[6] 



[7] 
[8] 

[9] 
[10] 

[11] 

[12] 
[13] 

[14] 

[15] 

[16] 

[17] 
[18] 
[19] 
[20] 
[21] 



[22] 
[23] 



[24] 
[25] 

[26] 



Present address: The Rudolf Peierls Centre for Theoret- 
ical Physics, 1 Keble Road, Oxford, 0X1 3NP, UK 
A.J. Liu and S.R. Nagel, Nature 396, 21 (1998). 
H.M. Jaeger, Physics World 18 (12), 34 (2005). 
C.S. O'Hern et al, Phys. Rev. E 68, 011306 (2003). 
C.S. O'Hern et al, Phys. Rev. Lett. 86, 111 (2001). 
T.G. Mason, J. Bibette and D.A. Weitz, Phys. Rev. Lett. 
75, 2051 (1995). 

For frictional packings, the situation is more complicated: 
frictional packings are not automatically marginal solids 
at jamming, but a critical-like scaling does emerge in the 
large-friction limit, see [l9l |. 
S. Alexander, Phys. Rep. 296, 65 (1998). 
M. Wyart, S.R. Nagel, and T.A. Witten, Europhys. Lett. 
72, 486-492 (2005). 

M. Wyart et al, Phys. Rev. E 72, 051306 (2005). 
M. Wyart, Ann. Phys. Fr. 30, 1-96 (2005). 
J. Geng et al, Phys. Rev. Lett. 87, 035506 (2001); Phys- 
ica D 182, 274 (2003); 

D. Serero et al, Eur. Phys. J. E 6, 169 (2001). 

A. Kasahara and H. Nakanishi, Phys. Rev. E 70, 051309 

(2004). 

C. Goldenberg and L Goldhirsch, Nature 435, 188 

(2005). 

N. Gland, P. Wang and H.A. Makse, Eur. Phys. J. E 20, 

179 (2006). 

I. Goldhirsch and C. Goldenberg, Eur. Phys. J. E 9, 245 

(2002). 

E. Somfai et al, Phys. Rev. E 72, 021301 (2005). 
A. Tanguy et al, Phys. Rev. B 66, 174205 (2002). 

E. Somfai et al, cond- mat/0510506 

F. Leonforte et al, Phys. Rev. B 70, 014203 (2004). 
W.G. EUenbroek et al, in Powders and Grains 2005, 
edited by R. Garcia-Rojo et al (A. A. Balkema, Rotter- 
dam, 2005), p. 377; and to be published. 

H.A. Makse et al, Phys. Rev. Lett. 83, 5070 (1999); 
Phys. Rev. E 70, 061302 (2004). 

For a floppy mode distortion, the second term in (1) can- 
cels the linear term —Yl,ij fij''^\\,ij ^^ ^^^ energy expan- 
sion. In (1) this term has been left out on account of force 
balance in the starting state. 

L.E. Silbert et al, Phys. Rev. E 65, 031304 (2002). 
L.E. Silbert, A.J. Liu and S.R. Nagel, Phys. Rev. E 73, 
041304 (2006). 

A very small fraction (typically 1%) of the responses is 
so disordered that the overlap with continuum behavior 
is minute (often because a cluster of particles far away 
shifts dramatically) — these have been omitted. 



